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FOREWORD 


This  report  documents  a  computer 
beamforming  process  of  a  planar  array. 


program  which  simulates 


The  program  uses  the  physical  characteristics,  electrical 
shading  characteristics,  medium  sound  speed  and  vehicle  motion  to 
obtain  detailed  Information  of  beam  structures.  Amplitude  and  phase 
coefficients  can  be  varied  to  study  main  and  side  lobe  structures  for 
any  combination  of  steered  or  unsteered  and  spread  or  unspread  beam 
patterns. 

i 

The  computer  program  documented  in  this  report  represents  the 
first  phase  of  a  beamforming  package,  and  Is  used  to  determine  the 
theoretical  limits  of  a  beamformer.  Work  Is  currently  under  way  on 
the  second  phase  of  this  study,  where  the  effects  of  parameter 
variation  during  production  will  be  examined.  A  third  phase  of  this 
study  will  deal  with  determining  beamformer  performance  In  a 
counter-measure  environment. 
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INTRODUCTION 


The  computer  program  THRVPLT  is  a  simulation  of  the  response  of 
a  desired  transducer  array.  It  is  currently  compatible  with  the  DEC 
VAX-11/780  computer-,  using  VAX-11  FORTRAN  code.  The  program  may  be 
used  to  generate  data  for  both  transmit  and  receive  beam  patterns. 
Calculations  are  made  relative  to  spherical  coordinates,  using 
necessary  coordinate  transformations  as  required.  Amplitude  (in 
decibels)  and  phase  angle  (In  degrees)  data  are  generated  for  a  beam 
over  any  desired  section  of  a  sphere.  The  input  data  created  for 
this  report  and  sample  outputs  are  Included  in  Appendices  A  and  B, 
respectively. 


£ 

iV- 


In  developing  the  program-,  several  assumptions  were  made.  The 
module  directivity  was  assumed  to  be  Independent  of  frequency,  depth 
and  location  In  the  array.  The  program  Included  in  this  report  is 
currently  set  up  for  planar  arrays,  and  uses  the  assumption  that 
there  will  be  no  variance  In  the  "Z"  direction  of  the  array  (see 
Figure  1).  That  Is,  the  face  of  the  array  will  be  smooth  across  its 
entirety.  Mutual  coupling  for  the  array  Is  set  to  zero.  Finally-, 
voltage  levels  below  -50  decibels  (db)  are  considered  to  be  noise, 
and  a  lower  boundary  is  set  at  that  level. 


The  program  has  a  set  of  run  options  available.  The  user  may 
specify  that  his  data  be  calculated  relative  to  the  main  response 
axis  (MRA),  to  the  "Z"  axis,  or  to  any  other  axis  he  chooses  (see 
Figure  1).  Once  this  has  been  chosen,  the  user  may  obtain  data  for  a 
variety  of  regions.  The  user  may  specify  that  his  data  be  calculated 
for  a  specific  roll  or  pitch  plane,  or  Input  a  range  of  pitch  and 
azimuthal  angles.  This  Is  discussed  In  detail  later.  With  respect 
to  module  directivity-,  the  user  must  Insert  Into  the  program  the 
directivity  pattern  he  chooses,  or  may  set  a  flag  to  run  an 
omnidirectional  pattern.  The  response  levels  may  be  normalized  to 
the  value  of  the  pattern  at  the  MRA,  to  the  sum  of  the  amplitude 
shading  coefficients-,  or  to  any  other  value  the  user  feels  Is  needed. 


Calculations  are  Included  In  the  program  for  equivalent  beam 
patterns  and  the  directivity  Index.  The  equivalent  beam  patterns  are 
calculated  assuming  Identical  transmit  and  receive  beam  patterns,  and 
should  be  viewed  as  a  special  case. 


*  w.  v 
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THEORY 

The  coordinate  system  used  In  this  development  Is  a  left-handed 
spherical  system,  as  shown  In  Figure  1.  This  system  was  chosen  to 
allow  calculations  to  be  made  relative  to  the  standard  azimuth  and 
pitch  angles.  Note  Is  made  that  the  angles  theta  (0)  and  phi  (0)  are 
not  the  usual  spherical  coordinates  as  found  In  most  text  books. 

The  major  calculations  made  In  the  program  yield  a  voltage 
response  (which  Is  later  converted  to  decibels)  and  a  phase  angle. 
The  program  assumes  that  waves  Incident  on  the  array  are  planar. 
This  assumption  Is  valid  providing  that  the  array-to-target  distance 
Is  much  greater  than  the  module  spaclngs  In  the  array. 

* 

A  plane  wave  may  be  represented,  using  complex  notation,  as 
*  *  A*EXP[1*(wt+B(l  ,m)+T(l-,m))] 


where 


B ( 1 ,m) 


amplitude 
2*ir*f  requency 
time 

phase  of  the  wave  at  the  (1  ,m)  element  relative  to  the 
array  center 


T ( 1  ,m)  *  electronic  phase  delay  of  the  (l,m)  element 


The  use  of  "EXP(x)"  In  the  equations  In  this  report  Is 
equivalent  to  the  natural  exponential  function  of  x.  The  term  "1"  In 
the  equations  Is  equal  to  the  square  root  of  negative  one. 

Using  a  maximum  amplitude  of  one*  and  letting  the  Individual 
module  responses  lie  between  zero  and  one,  the  response  of  the  (1,m) 
element  may  be  written  as 

V ( 1 , m )  »  H(1 ,m)*R(l ,m)*EXP[1*(wt+B(l ,m)+T(l ,m))] 

where  R(1,m)  and  H(1,m)  are  the  electrical  attenuation  and  module 
directivity-,  respectively',  of  the  (1  ,m)  element  (both  are  normalized 
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to  one).  A  discussion  of  module  directivity  Is  Included 
this  report.  Neglecting  the  time  dependence  of  the 
response  may  be  written  as 

V{ 1 ,m)  =  H{1  ,m)*R(l ,m)*EXP[i*(B( 1 ,m)+T(l ,m))] 


later  In 
wave,  the 


The  first  phase  term,  B(1,m),  is  simply  the  difference  in  phase 
at  which  the  wave  arrives  at  element  (1 ,m)  with  respect  to  the  array 
center.  If  the  phase  term  B(l,m)  has  a  positive  value,  it  is  said  to 
be  a  phase  lead,  while  a  negative  value  of  B(l,m)  Is  defined  as  a 
phase  lag.  This  phase  delay  relates  directly  to  the  perpendicular 
distance  between  the  plane  wave  front  and  the  element  being  looked  at 
when  the  wave  front  Is  at  the  array  center.  This  distance  (phase 
delay)  can  be  seen  by  examining  Figures  2  and  3. 


Figure  2  is  a  sketch  of  a  plane  wave  front  Incident  on  an  array, 
perpendicular  to  the  direction  given  by  theta  and  phi.  The  wave 
front  Is  represented  by  the  broken  line.  Module  1  Is  located  at  the 
center  of  the  unprimed  axes,  while  module  2  is  located  at  the  center 
of  the  primed  axes.  The  origin  (0,0,0)  of  the  unprimed  axes  will  be 
taken  as  the  reference  point  for  the  derivation  that  follows.  The 
line  segments  from  module  1  to  point  A,  point  A  to  point  B,  and  point 
B  to  module  2  represent  the  distances  X(k),  Y  ( k )  and  Z(k), 
respectively.  The  origin  of  the  primed  axes  Is  located  at  the  point 
(X(k ) , Y(k)  ,Z(k ) ) .  The  vector  S  represents  the  total  perpendicular 
distance  to  the  wave  front  from  module  2. 

To  find  the  total  distance  to  the  wave  from  module  2,  the 
modules  will  first  be  looked  at  as  If  they  were  both  In  the  X-Z 
plane,  with  module  1  being  located  at  (0,0,0)  and  module  2  being 
located  at  (X(k) ,0,Z(k) ).  The  contribution  due  to  a  difference  in 
the  Y  positions  will  then  be  added. 

To  start',  let  the  plane  wave  be  Incident  on  the  modules  in  the 
direction  (9,0),  l.e.-,  where  9*9  and  0=0.  This  Is  shown  In  Figure 
3(a).  As  shown  in  Figure  3(a),  the  value  of  X(k)  Is  positive,  while 
Z ( k )  Is  negative.  The  value  of  theta  as  shown  Is  also  negative. 
With  these  conditions,  the  total  distance  to  the  wave  front  from 
module  2  Is 

51  *  X(k)*s1n(9)  +  Z (k )*cos(9) 

Moving  to  figure  3(b),  there  Is  still  no  Y  separation,  but  the 
wave  Is  now  Incident  In  the  direction  (9,9).  The  view  In  Figure  3(b) 
Is  taken  looking  down  the  wave  front  at  a  vertical  plane.  With  the 
wave  now  Incident  In  the  direction  (9,0),  the  total  distance  to  the 
wave  front  from  module  2  becomes 

52  -  Sl*cos(0)  -  [X(k)*s1n(9)  +  Z(k)*cos(9)]*cos(0) 
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Figure  3(c)  shows  the  representati on  of  the  addition  of  a  Y 
separation  for  the  modules.  The  view  is  the  same  as  that  used  in 
Figure  3(b).  In  this  figure,  module  2  has  been  moved  a  distance  Y ( k ) 
below  module  1  (l.e.,  Y(k)  has  a  negative  value)',  and  the  wave  is 
again  Incident  In  the  direction  (6,0).  Modules  1  and  2  are  now 
located  at  the  positions  (0,0,0)  and  (X(k)  ,Y(k) , Z ( k ) ) ,  respectively. 
The  total  distance  to  the  wave  front  from  module  2  is  then 


S  -  X(k)*s1n(0)*cos(0)  +  Y(k)*s1n(0)  +  Z (k )*cos(9)*cos(0) 

In  the  equation  representing  V(l,m),  B(l,m)  is  In  fact  directly 
related  to  this  distance  S.  Adding  a  normalizing  term  (to  convert  S 
to  a  phase)  and  letting  the  subscript  "k"  become  "l,m",  B(l,m)  may  be 
written  as 


B(l,m)  -  (2VX)*[X(1  ,m)*s1n(6)*cos(0) 

+  Y  ( 1 ,m)*s1n(0)  +  Z(l  ,m)*cos(0)*cos(0)3 
•  (2*f/c)*[X(l ,m)*s1n(0)*cos(0) 

+  Y ( 1  ,m)*s1n(0)+Z>(1  ,m)*cos(9)*cos(0)] 

where 


X(l,m) 

-  X 

location 

of 

(1  ,m) 

el ement 

rel atl ve 

to 

(0,0,0) 

Y(l,m) 

«  Y 

location 

of 

(1  ,m) 

element 

relative 

to 

(0,0,0) 

Z(l,m) 

-  Z 

location 

of 

(1  ,m) 

el ement 

rel atl ve 

to 

(0,0,0) 

X  »  wavelength 

c  ■  speed  of  wave  In  medium 
f  ■  frequency 


The  second  phase  term',  T(l,m)',  Is  the  steer/spread  coefficient. 
The  program  will  calculate  these  coefficients  for  a  steered  beam,  but 
not  for  a  spread  beam,  using  the  equation 

T ( 1 , m )  -  -W*[X(1 ,m)*s1n(0(MRA))*cos(0(MRA)) 

+  Y(1 ,m)*s1 n(0(MRA) )+Z ( 1 ,m)*cos (0(MRA ) ) *cos (0(MRA ) ) ] 

where  (0(MRA))  and  ( 0 ( MR A ) )  are  the  values  of  theta  and  phi  at  the 
MRA.  The  term  W  above  Is  equal  to  2irf/c,  the  normalization  term. 
Phase  coefficients  containing  spread  Information  must  be  calculated 
by  the  user  before  running  the  program.  The  Z(l,m)  terms  have  been 
omitted  In  this  simulation,  using  the  assumption  that  the  face  of  a 
planar  array  will  be  smooth  across  Its  entirety,  if  the  user  foresees 
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a  Z  dependence  In  his  system,  it  should  ;  >ut  into  the  equations  in 
the  program.  The  frequency  used  in  these  equations  may  be  adjusted 
for  motion  using  the  common  Doppler  shift  equation,  as  found  in  texts 
on  general  physics.  The  equation  Is  included  in  the  program  and  may 
be  flagged  as  desired. 

Using  Euler's  theorem,  and  letting  e(l ,m)  *  T ( 1 ,m)+B(l  ,m) ,  the 
response  may  be  written  as 

V(l,m)  =  H(1 ,m)*R(1 ,m)*EXP[i*(e(l ,m))] 

*  H (1  ,m)*R (1  ,m)*[cos(e(l ,m))  +  i*sin(e(i ,m))] 

To  get  the  total  response  of  the  array,  the  individual 
transducer  responses  are  summed  up.  This  total  response  is 

Ny  Nx 

Vs*  E  E  H(1  ,m)*R(l ,m)*[cos(e(l ,m) )+i*sin(e(l ,m) )] 

£-1  m-l 

where  Nx  and  Ny  are  the  number  of  columns  and  rows,  respectively. 

Taking  the  modulus  of  the  complex  term  above,  the  response 
becomes 

»y  Nx  , 

V*([Z  E  H(l,m)*R(l,m)*cos(e(l,m))]2 
l-l  m-l 


Ny  Nx 


+  [  E  E  H(l,m)*R(l,m)*s1n(e(l,m))r) 

t-1  m-i 


2  .  1/2 


The  assumption  Is  now  made  that  the  term  H(l,m)  will  be  the  same  for 
every  element;  this  Is  discussed  later  In  this  section.  Using  this 
assumption,  and  letting  all  H(l,m)  be  equal  to  H,  the  response  may  be 
written  as 

Ny  Nx  , 

v  -  H*( [  E  E  R(l,m)*cos(e(l,m))3z 

l-l  m-l 


Ny  Nx  , 

+C  E  E  R(l  ,m)*s1n(e(l  ,m) )]  z) 

l-l  m-l 


2.1/2 


The  phase  of  the  total  response.  In  the  complex  system.  Is  given 
by  the  vector  sum  of  the  Individual  phases.  Specifically,  the  phase 
Is  equal  to 

e  ■  arctanClm/Re"1 


/  A 
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where 


*  phase  of  the  total  response 
=  sum  of  the  Imaginary  terms 

*  sum  of  the  real  terms 


The  assumption  was  made  earlier  that  the  module  directivity  will 
be  the  same  for  all  elements.  This  assumption  has  been  shown  through 
hardware  testing  to  be  a  reasonable  one  for  the  systems  looked  at 
thus  far.  For  values  of  theta  and  phi  out  to  +60  degrees, 
differences  in  directivity  between  elements  were  negligible.  Caution 
is  given  here  that  this  assumption  may  not  be  valid  for  other  arrays 
of  Interest.  It  Is  left  to  the  user  to  test  his  system  to  see  if  the 
assumption  holds',  and  to  make  program  modifications  if  necessary. 
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PROGRAM  DESCRIPTION 


The  main  section  of  the  program  is  controlled  by  two  nested 
loops.  The  limits  for  these  loops  are  based  on  the  spherical 
coordinates  theta  (0)  and  phi  (0)-,  as  shown  in  Figure  1.  The  inner 
loop  runs  over  the  desired  range  of  theta,  while  the  outer  loop  runs 
over  the  desired  range  of  phi.  Both  loops  must  be  run  with  a 
constant  increment. 

The  user  may  run  the  program  relative  to  any  position  he 
chooses.  Usual  positions  to  run  relative  to  are  the  Z  axis 
(theta  =  phi  =  0)  or  the  main  response  axis  (MRA ) . ’  When  running 
relative  to  an  axis  other  than  the  Z  axis,  adjustments  are  made  to 
theta  and  phi,  through  the  use  of  offsets. 

When  running  the  program,  the  user  must  specify  the  set  of 
angles  he  is  interested  in.  He  may  do  this  directly,  or  may  specify 
a  roll  plane  and  range  of  angles  relative  to  the  roll  axis.  The  roll 
axis  is  the  ray  about  which  the  plane  is  rotated.  The  roll  angle  is 
to  be  given  between  -90  and  90  degrees  (e.g.,  specify  -70  instead  of 
110).  In  Figure  C-l  (in  Appendix  C),  the  Z  axis  is  the  roll  axis. 
When  the  user  specifies  a  roll  plane  to  be  calculated,  the  program 
converts  back  to  real-space  theta  and  phi.  The  angle  on  the  roll 
plane,  relative  to  the  roll  axis,  will  be  designated  here  as  theta 
prime  (O'),  while  the  roll  angle  itself  will  be  designated  gamma  (Y). 
The  following  relations  are  used  for  the  conversion  from  the  roll 
plane  to  real-space  coordinates: 

9  =  arctan[tan(9* )  *  cos(Y)]  +  theta  offset 

9  9  arcsin[s1n(9' )  *  sin(Y)]  +  phi  offset 

A  derivation  and  diagram  of  these  relations  is  included  in  Appendix  C 
(see  Figure  C-l). 

The  total  response  is  normalized  to  a  value  requested  by  the 
user.  While  any  value  may  be  chosen,  the  most  common  and  useful 
choices  are  either  to  the  total  system  or  to  the  individual  beam 
Itself  (l.e.,  to  the  response  at  the  MRA).  System  and  self¬ 
normalizations  can  be  calculated  in  the  program  by  setting  the  flag 
appropriately;  any  other  choice  for  normalization  must  be  input  by 
the  user. 
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When  the  total  normalized  response  and  phase  have  been 
calculated,  minor  adjustments  are  made.  The  response  is  converted  to 
a  decibel  level,  with  a  lower  limit  set  to  -50  db.  The  phase  angle 
is  calculated  using  a  four-quadrant  arctangent  function,  and  will 
have  a  value  between  +180  degrees. 

As  mentioned  earlier,  the  program  will  take  into  account  shaping 
of  the  beam  due  to  module  directivity.  The  equation  for  this  shaping 
must  be  put  into  the  program  by  the  user-,  and  should  be  normalized  to 
one.  The  variable  representing  this  factor  is  designated  "H",  and  is 
found  in  both  the  normalization  and  main  sections  of  the  program. 

The  shaping  factor  is  employed  as  a  multiplier  on  the  voltage 
response  levels.  The  shaping  factor  currently  in  the  model  is 

H=[.6+cos(9)cos(0)] 2/2. 56 

The  equation  above  is  an  approximation  to  a  cardioid  pattern  and 
will  have  a  value  between  zero  and  one.  Figure  4  shows  a  polar  plot 
of  this  pattern.  This  directivity  pattern  is  symmetric  about  the  Z 
axis,  and  is  shown  in  Figure  4  in  the  X-Z  plane.  It  is  not  to  be 
considered  accurate  in  the  region  directly  behind  the  face  of  the 
array.  When  running  the  program,  a  flag  may  be  set  to  default  to  an 
omnidirectional  pattern  (i.e.-,  H=1  over  the  entire  space). 

Additional  calculations  are  performed  in  the  program  to  give 
levels  for  the  directivity  index  and  equivalent  beam  patterns  for  use 
in  reverberation  calculations.  The  equivalent  beam  patterns  here  are 
calculated  assuming  Identical  transmit  and  receive  beams.  This 
should  be  taken  Into  account  when  examining  the  output. 
Additionally,  these  values  will  only  be  valid  when  the  program  is  run 
over  the  entire  necessary  space.  For  example,  to  get  a  good  value 
for  the  directivity  index,  the  total  sphere  should  be  run.  Running 
over  the  entire  sphere,  however,  requires  knowledge  of  the  beam 
pattern  and  module  directivity  for  that  entire  region.  A 
hemispherical  run  using  0  and  0  between  +90  degrees  has  in  general 
proven  to  be  satisfactory  for  a  planar  array,  since  contributions 
from  behind  the  array  are  usually  quite  small.  The  angular 
resolution  Is  left  up  to  the  user. 

The  following  Integrals  are  Included  in  the  program: 

0(2)  <p(2) 

D I  =  10*1  og (4tt )  -  10*1  og[  /  /  b(0,0)*cos(0)d0d0] 

6(1)  <J>(1) 

0(2)  <J>(2) 

EBPV0L  «  10*1 og[  /  /  b(0,0)*b' (0,0)*cos(0)d0d0] 

0(1)  4>(1) 

0(2) 

EBPSB  -  10*1 og[  /  b(0 ,0 (MRA ) ) *b ' (0 ,0 (MRA )  )d0] 
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where  DI  is  the  directivity  index  and  EBPVOL  and  EBPSB  are  the  volume 
and  surface/bottom  equivalent  beam  patterns,  respectively.  The 
transmit  and  receive  patterns  are  given  by  b'(0,0)  and  b(0,0)» 
respectively.  The  program  currently  automatically  sets  the  value  of 
phi  for  the  surface  integration  to  the  value  of  phi  at  the  MRA.  The 
limits  on  the  integrals  are  set  by  the  user.  Theoretically,  the 
limits  should  cover  the  entire  sphere.  However,  as  mentioned 
earlier,  a  hemisphere  has  generally  been  found  to  be  satisfactory. 

The  integrations  in  the  program  are  done  numerically',  using  a 
rectangular  integration  method.  A  discussion  of  the  method  is 
included  in  Appendix  D.  As  an  example  of  the  method's  accuracy,  the 
directivity  index  was  calculated  by  the  program  for  an 
omnidirectional  voltage  response  and  a  cos(0)*cos(0)  voltage  reponse. 
Direct  integrations  were  also  done  for  these  same  integrals,  and  the 
resultant  errors  were  calculated.  The  results  are  given  below  in 
Table  1. 


TABLE  1.  COMPARISON  OF  DIRECTIVITY  INDEX  (DI)  CALCULATIONS 


Omnidi rectional 

Computer  Calculation  0.0003  db 

Direct  Integration  0.0000  db 


cos(0)*cos(fil) 
4.7713  db 
4.7712  db 


The  error  in  the  omnidirectional  response  is  approximately 
.0003  db',  while  the  cos(0)*cos(jJ)  response  has  an  error  of 
approximately  .0001  db.  The  integrations  on  the  computer  were  done 
using  a  resolution  of  one  degree  in  each  direction.  The  error  may  be 
attributed  to  resolution  error,  and  for  most  applications  is 
negl 1 g 1 b 1 e . 

Table  2  below  contains  four  sample  calculations  of  equivalent 
beam  patterns.  All  four  patterns  were  calculated  using  the  same 
receive  beam  (b(0,0)  -  see  Appendix  A  for  inputs)',  while  the  transmit 
pattern  (b'(0,(3))  was  varied.  For  cases  3  and  4,  the  transmit 
patterns  are  cone  shaped;  within  the  cone  specified,  the  pattern  has 
the  form  of  a  zero  decibel  omnidirectional  beam,  while  outside  the 
cone  it  is  assumed  to  have  a  negligible  level. 


TABLE  2.  SAMPLE  EQUIVALENT  BEAM  PATTERN  (EBP)  CALCULATIONS 


1.  b'(0,0)  *  b(0,0) 

a.  EBPVOL  =  -14.7485  db 

b.  EBPSB  =  -7.3586  db 


9 


2.  b'(9,0)  =  zero  db  omnidirectional 

a.  EBPVOL  »  -11.7097  db 

b.  EBPSB  =  -5.8553  db 

3.  b * (9,0)  ■  zero  db  omnidirectional  out  to  receive 

pattern's  -10  db  points 

a.  EBPVOL  *  -12.1470  db 

b.  EBPSB  *  -6.0466  db 


4.  b'(9,j?)  *  zero  db  omnidirectional  out  to  receive 

pattern's  -3  db  points 

a.  EBPVOL  =  -13.6843  db 

b.  EBPSB  *  -6.8339  db 


Case  3  above  approximates  a  constant  transmit  source  level 
within  the  10  db  receive  band  width  with  the  transmit  beam  falling 
off  below  -50  db  outside.  The  result  given  In  case  3  above  is  a  good 
representation  of  the  real-world  system,  since  transmit  beams  will 
usually  be  shaped  somewhere  between  a  shaded  beam  and  an 
omnidirectional  beam.  By  comparing  cases  1  through  3,  It  Is 
reasonable  to  assume  that  when  the  source  level  Is  constant  over  the 
10  db  receive  band  width,  the  real-world  equivalent  volume  beam 
pattern  will  be  2.6  db  to  3  db  higher  than  the  program's  calculation, 
while*  the  equivalent  surface/bottom  beam  pattern  will  be  1.3  db  to 
1.5  db  higher  than  the  program's  calculation.  Case  4  points  out  that 
us1n?  the  receive  pattern's  -3  db  points  as  a  guide  will  give  rise  to 
an  appreciable  difference-,  and-,  depending  on  requirements-,  may  or  may 
not  be  an  acceptable  approximation. 
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SUMMARY 


The  computer  program  documented  In  this  report  is  useful  for 
determining  the  response  of  planar  arrays.  It  may  be  used  as  an  aid 
In  determining  the  combination  of  amplitude  and  phase  shading 
coefficients  which  will  produce  a  desired  beam  pattern.  The  program 
will  also  prove  useful  in  determining  the  effects  of  variations  in 
frequency  and  wave  propagation  speed. 

The  results  from  the  program  will  represent  the  theoretical 
limits  of  an  array.  Actual  performance  will  most  likely  be  somewhat 
degraded  in  production  units,  particularly  in  side  lobe  structures. 


11 


NSWC  TR  82-376 


FIGURE  2.  INCIDENT  PLANE  WAVE  IN  THREE  DIMENSIONS 


v  _i,  _v - 
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FIGURE  3(a).  TWO  MODULES  LOCATED  AT  (0,0,0)  AND  (X(k), 
PLANE  WAVE  INCIDENT  IN  DIRECTION  OF  (0,0) 


FIGURE  3(b).  TWO  MODULES  LOCATED  AT  (0,0,0)  AND  (X(k), 
PLANE  WAVE  INCIDENT  IN  DIRECTION  OF  (0,0) 


,Z(k)) 


,Z(k)) 
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FIGURE  3(c).  TWO  MODULES  LOCATED  AT  (0,0,0)  AND  (X(k ) ,  V(k  )  ,2 (k  ) ) 
PLANE  WAVE  INCIDENT  IN  DIRECTION  OF  (0,0) 
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APPENDIX  A 
INPUT  VARIABLES 


Inputs  to  the  program  are  to  be  put  into  a  data  file  called 
THRYPLT.  The  program  references  this  file  for  all  of  its  inputs, 
excluding  matrices.  The  program  uses  four  matrices  which  must  be 
created  by  the  user:  X  and  Y  positions',  amplitude  shading  and  phase 
shading.  The  X  and  Y  position  matrices  give  module  locations.  In 
Inches*,  relative  to  the  array  center.  The  amplitude  shading 
coefficients  are  given  relative  to  a  maximum  value  of  one,  and  the 
phase  shading  coefficients  are  given  In  degrees.  A  Z  position  matrix 
will  be  requlred-if  the  user  foresees  any  Z  dependence  in  the  array. 
The  Input  data  file  and  matrices  used  to  generate  this  report's 
output  data  are  given  In  Figure  A-l.  Below  is  a  list  of  the  Input 
variables  and  their  definitions.  Numbers  to  the  left  of  the  variable 
names  correspond  to  the  lines  of  data  under  INPUT  DATA  FILE  in  Figure 
A-l.  Input  numbers  2,  3*,  4  and  6  are  Input  with  a  character  format. 
All  of  the  other  Inputs  use  a  free  format,  and  decimal  points  are  not 
required  for  real  variables. 


1.  Phaseplot 


2.  FI  1  x 

3.  FI  ly 

4.  FI  1  r 

5.  Tscalc 


6.  Flit 


7.  Loop  Control 

a.  Lai 

b.  La2 

c.  Inca 

d.  Lbl 

e.  Lb2 

f.  Incb 

g.  Aresol 

h.  Bresbl 


Controls  which  plots  will  be  made 
0  ■  db  vs.  angle 

1  *  db  vs.  angle  and  phase  vs.  angle 

2  ■  phase  vs.  angle 

3  ■  no  plots  made 

Name  of  the  X  position  matrix 
Name  of  the  Y  position  matrix 
Name  of  the  amplitude  shading  matrix 
Controls  If  phase  matrix  Is  to  be 
Input  or  calculated  by  the  program 
0  -  Input  a  matrix  name 
1  ■  program  will  calculate  steering 
coefficients 

Name  of  the  phase  shading  matrix  - 
present  only  If  #5  Is  a  0 

Input  a  1 

Number  of  phi  angles  to  be  run 
Input  a  1 
Input  a  1 

Number  of  theta  angles  to  be  run 
Input  a  1 

Resolution  for  phi  loop  (In  degrees) 
Resolution  for  theta  loop  (In  degrees) 
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i.  As tart 

j.  Bstart 


Starting  position  for  phi  (in  degrees) 
Starting  position  for  theta  (in  degrees) 


Note:  Some  loop  control  variables  are  used  in  the  output 
sections  of  the  program.  The  default  values  above 
should  be  used  to  insure  correct  output  formats. 


8.  a.  F 

b.  C 

c.  Nx 

d.  Ny 

9.  a.  Vweap 

b.  Amalnr 

c.  Bmainr 

d.  Norml 


e.  Freqfl 


10.  Anorml 

11.  Offlag 


12.  Frol  1 


13.  Roll 

14.  Hflag 


15.  FD.CD 


Frequency  of  the  wave  (in  kHz) 

Speed  of  wave  in  medium  (in  ft/sec) 

Number  of  columns  in  the  array 
Number  of  rows  in  the  array 
Velocity  of  the  weapon  (in  knots) 

Value  of  phi  (in  degrees)  at  the  MRA 
Value  of  theta  (in  degrees)  at  the  MRA 
Normalization  flag 

0  =  normalization  is  equal  to  the 
sum  of  the  amplitude  shading 
coefficients 

1  =  normalize  to  the  beam  Itself 

2  =  normalize  to  an  input  value 

Controls  if  angular  doppler  shift  is 
to  be  added  Into  the  frequency 

0  ■  no  ' 

1  ■  yes 

Normalization  value  (In  decibels)  > 
present  only  if  #9-d  Is  a  2 
Controls  if  beam  is  to  be  offset  to  the  MRA 
0  *  no 

1  *  yes 

Controls  If  a  roll  plane  calculation 
Is  requested 
0  =  no 

1  *  yes 

Roll  plane  angle  (In  degrees)  - 

present  only  If  #12  is  a  1 

Controls  the  type  of  directivity  shaping 

present 

0  ■  the  equation  In  the  program  Is  used 
1  ■  an  omnidirectional  pattern  Is  used 
Design  frequency  (kHz)  and  wave 
speed  (ft/sec).  Need  to  be  present 
only  If  #5  Is  a  1. 
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APPENDIX  B 
OUTPUT  FILES 


The  user  has  the  option  of  selecting  his  outputs  from  the 
following  list: 

1.  A  formatted  file  combining  an  input  list  and  selective 
outputs  such  as  the  directivity  index. 

This  file  is  found  under  F0R003.DAT  (see  Figure  B-l). 

2.  A  sequential  file  of  decibel  and  phase  data,  sequenced 
according  to  the  phi  and  theta  loops. 

This  file  is  found  under  F0R013.DAT. 

3.  A  formatted  file  divided  into  two  sections: 

a.  decibel  values  vs.  angle  (see  Figure  B-2) 

b.  phase  angles  vs.  angle  (see  Figure  B-3) 

These  are  found  under  F0R002.DAT. 

4.  A  file  containing  plots  of  the  output  data: 

a.  decibel  values  vs.  angle  (see  Figure  B-4) 

b.  phase  angles  vs.  angle  (see  Figure  B-5) 

These  plots  are  found  under  PLOTDAT.DAT 


The  angles  used  in  the  tables  of  output  number  3  are  for  a 
specific  plane  requested  by  the  user.  For  example,  if  the  user 
specifies  a  pitch  plane,  the  angles  represent  the  values  of  0.  If  a 
roll  plane  Is  requested,  the  angles  in  the  tables  represent  values  of 
0'.  Output  number  3  should  be  used  in  conjunction  with  number  4  if 
the  accuracy  of  the  plots  is  insufficient.  When  no  plots  are 
requested  (l.e.-,  input  number  one  is  set  to  3),  the  FOR002.DAT  file 
has  a  slightly  different  format  than  above.  It  will  contain  all  of 
the  decibel  and  sign  data  generated  by  the  program,  and  will  have  the 
pitch  angles  labeled  at  the  top.  However,  this  file  will  contain 
only  the  decibel  level  and  sign  of  the  real  part  of  the  voltage 
response.  It  will  not  contain  the  actual  phase  angle  data  in  this 
mode  (see  Figure  B-6).  When  using  the  program  to  generate  plots,  the 
F0R002.DAT  file  will  contain  the  last  set  of  decibel  and  phase  data, 
as  seen  in  Figures  B-2  and  B-3. 

To  use  the  plot  routines  In  the  program,  the  user  must  have 
access  to  a  Digital  DECwriter  III  or  similar  hardcopy  terminal.  In 
addition,  the  user  must  create  a  data  file  named  PLOTDAT,  having  an 
available  record  length  of  at  least  200.  This  is  to  be  done  using  a 
file  create  routine  from  the  computer,  or  In  the  OPEN  statement  for 
that  file  In  the  program.  The  attached  program  code  currently  uses  a 
file  created  by  a  file  create  routine.  To  print  the  plots  out',  the 


B-l 


following  settings  must  be  made  on  the  terminal: 

1.  Horizontal  print  control  is  set  to  16  characters/inch. 

2.  Vertical  print  control  is  set  to  8  lines/inch. 

3.  Page  size  is  set  to  88  lines/page. 

The  required  form  size  is  approximately  15  inches  wide  by  11 
inches  long.  The  user  must  also  set  the  terminal  width  to  200.  On 
the  VAX  11/780,  the  command  which  will  set  this  is 

'SET  TERMINAL/WIDTH=200‘. 

The  above  settings  should  also  be  used  if  the  user  desires  that  the 
file  F0R002.DAT  fit  into  a  standard  8-1/2  by  11  inch  notebook. 
Reductions  are  necessary  for  the  plots  to  fit  into  these  same  books. 

When  using  the  data  plots,  note  should  be  taken  of  the  degree  of 
accuracy  inherent  in  a  computer;  i.e.,  there  is  a  finite  limit. 
When  using  the  phase  vs.  angle  plot,  it  is  seen  that  the  values  of 

the  phase  are  between  _+ 1 80  degrees.  In  terms  of  phase,  -180  and  180 

degrees  are  the  same.  However,  the  computer's  finite  limit  in 
accuracy  may  cause  the  phasor  to  oscillate  around  this  point,  and  be 
plotted  at  two  different  values. 

i 

When  using  the  sign  vs.  angle  plot  (located  under  the  db  vs. 

angle  plot),  the  "sign"  term  represents  the  sign  of  the  real  part  of 

the  total  voltage  response,  i.e.,  whether  It  is  positive  or  negative. 
This  tells  the  user  which  quadrant  set  (1-4  or  2-3)  the  phasor  is 
located  in  for  a  given  angle.  The  points  having  a  value  of  0  for 
sign  are  artificially  set  to  that  value  to  tell  the  user  that  the 
value  (in  decibels)  has  gone  below  -50  db. 

If  the  user  requests  the  program  to  calculate  the  phase  steering 
matrix,  it  will  be  output  to  a  file  named  F0R020.DAT.  The  user 
should  rename  this  file  to  reflect  its  contents.  For  example,  if  the 
file  contains  phase  steering  data  for  a  (10,10)  steer,  the  file  might 
be  renamed  to  T1010.DAT  to  remind  the  user  that  the  file  contains  "T" 
data  for  a  "(10,10)"  steer  angle. 
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FIGURE  B-l.  INPUT/OUTPUT  SUMMARY  FILE 
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FIGURE  B-2.  DECIBEL  VS.  ANGLE  FILE 


GECIKLS 


FIGURE  B-4.  DECIBEL  VS.  ANGLE  PLOT 


phase  mu  vs  ncTs 


FIGURE  B-5.  PHASE  VS.  ANGLE  PLOT 
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APPENDIX  C 

CONVERSION  FROM  ROLL  PLANE  TO  REAL-SPACE  COORDINATES 


Following  is  a  development  of  the  equations  used  to  convert  from 
a  given  roll  angle  and  theta  prime  to  the  real-space  coordinates 
theta  (0)  and  phi  (0).  The  angles  theta  and  phi  are  as  noted 
earl  ler. 

In  Figure  C-1-,  the  angle  gamma  (T)  is  the  roll  plane  angle  and 
theta  prime  ( 0 * )  is  the  angle  from  the  axis  on  that  roll  plane. 

Using  Figure  C-l,  the  following  relationships  may  be  written 
down  for  the  four  angles  mentioned  above: 

1)  tan(O)  *  S/R 

1 1 )  sln(O)  -  S/P 

111)  cos ( 9 )  *  R/P 

1 v)  tan(0)  *  W/P 

v)  s1n(0)  *  W/Q 

vl)  cos(0)  »  P/Q 
vl 1 )  tan(y  )  »  W/S 

vl 1 1 )  sln(y)  «  W/T 
lx)  cos ( y)  *  S/T 

x)  taniQ' )  ■  T/R 
x 1  )  sln(Q')  -  T/Q 
x 1 1 )  cos(9' )  ■  R/Q 

where  the  letters  P,  Q,  R,  S,  T  and  W  are  the  lengths  designated  In 
Figure  C-l. 

Combining  equations  (1),  (lx)  and  (x)-,  the  following  relation 
may  be  seen: 

cos(Y)*tan(9‘ )  ■  S/T*T/R  -  S/R  »  tan(9) 


which  leads  to  the  relation 
9  *  arctan (tan (9 ' )*cos{Y) ) 


Combining  equations  (v),(v111)  and  (xi),  the  following  relation 
may  be  seen: 


si n( Y)*s1n(9 1 )  -  W/T*T/Q  -  W/Q  -  s1n(0) 
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which  gives 

9  *  arcs1n(s1n(9' )*sin(Y)) 


These  equations  are  then  used  In  the  program  to  make  the 
necessary  transformations.  Theta  and  phi  offsets  are  added  to  the 
theta  and  phi  equations,  respectively,  to  change  the  location  of  the 
roll  axis. 


1 
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APPENDIX  D 


INTEGRATION  TECHNIQUE 


The  Integrations  done  In  the  program  are  done  in  spherical 
coordinates,  using  Integrals  of  the  form 

0(2)  <(,(2) 

1  =  /  /  F*cos(0)d0d0 

0(1)  4>(D 

where  F  represents  a  beam  pattern  (or  product  of  beam  patterns)  and 
0(1),  0(2),  0(1)  and  0(2)  are  the  Integration  limits. 

The  computer  program  performs  the  Integrations  numerically  using 
the  Newton-Cotes  method  of  degree  zero.  A  resolution  of  one  degree 
has  normally  been  used,  along  with  the  mid-point  rule.  Any  other 
resolution  may  be  chosen;  thus  far,  one  degree  has  been  found  to  be 
satl sf actory . 

The  degree  zero  Newton-Cotes  method  sets  the  entire  interval 
used  to  the  value  at  the  Interval's  mid-point.  The  accuracy  of  this 
method  depends  directly  on  the  second  derivative  (f'‘)  of  the  curve 
over  the  intervals.  This  may  also  be  thought  of  as  the  derivative  of 
the  slope  over  each  Interval.  If  the  slope  over  an  interval  Is 
constant  (f''=0),  the  error  for  that  Interval  will  approach  zero, 
while  a  large  change  in  the  slope  over  the  Interval  will  result  in  a 
large  error  for  that  Interval.  For  the  two  beams  analyzed  In  Table  1 
(see  page  9)',  the  slopes  were  mainly  slow  changing-,  and  the 
one-degree  resolution  gave  a  very  small  error. 

For  the  sonar  beam  created  for  this  report,  a  large  value  of  f 1 1 
was  found  at  the  null  points.  However,  the  levels  at  these  points 
are  so  depressed  that  the  resulting  relative  error  Is  very  small. 
For  this  beam  (see  the  Inputs  In  Appendix  A),  two  separate 
directivity  Index  Integrations  were  run:  the  first  was  over  the 

entire  hemisphere,  while  the  second  used  only  those  points  contained 
In  the  main  lobe  (0  and  0  between  +17  degrees).  The  result  for  the 
hemispherical  run  was  22.  7025  d"F,  while  the  main  lobe  Integration 
yielded  a  value  of  23.0245  db.  For  this  beam-,  then,  the  side  lobe 
contribution  was  less  than  .33  db.  While  this  Is  one  discrete  case, 
It  Is  seen  that  eliminating  all  of  the  null  points  and  side  lobes  did 
not  drastically  change  the  result.  The  user  is  cautioned  to  check 
his  beams  as  he  creates  them  to  see  that  his  values  for  f ' '  are  small 
over  his  range  of  Interest.  To  Increase  the  degree  of  accuracy,  the 
user  may  reduce  the  Interval  size.  However,  a  decrease  In  the 
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interval  size  will  result  in  an  increase  In  run  time.  For  example, 
decreasing  the  Interval  size  from  one  degree  to  one  half  of  a  degree 
In  each  direction  will  result  in  an  increase  in  run  time  of  almost 
four  times. 

Figure  D-l  Is  a  representation  of  the  interval  used  in  the 
integrations.  The  delta  area  is  centered  at  (9,0).  R  represents  the 
Intensity  of  the  beam  at  (9,0).  The  area  can  be  approximated  as  a 
rectangle,  whose  sides  are  the  segments  labeled  2  and  3.  For  small 
values  of  d9  and  d0,  arc  lengths  can  be  approximated  by  their 
corresponding  chord  lengths.  Then  as  shown  in  Figure  D-l,  segment  3 
Is  equal  to  d0,  while  segment  2  is  approximately  equal  to  segment  1 
times  the  cosine  of  0',  namely  cos(0)d9.  Multiplying  segments  2  and  3 
gives  a  value  for  the  delta  area  of  cos(0)d0d9.  For  small  values  of 
d0',  surface  2  will  vary  less  than  .01  percent  from  (0-(d0/2))  to 
(0+(d0/2)),  making  the  approximation  used  quite  reasonable. 


D-2 
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APPENDIX  E 
PROGRAM  LISTING 


C  PROGRAM  8PATRN(INPUT, OUTPUT, TAPE2, TAPE3, TAPE9, TAPE10, TAPE11  , 

C  TAPE5.TAPE12) 

DIMENSION  X(8,8),T(8,8) , Y ( 8 , 8 ) ,R(8,8) 

DIMENSION  A ( 365 ) ,B(365) 

DIMENSION  VI (365) ,V2(  365) ,V3(365)  ,V4 ( 365 )  ,V5(365)  , 

*SS 1(365) ,SS2( 365) ,SS 3(365 ) ,SS4( 365)  ,SS5( 365) 

DOUBLE  PRECISION  TB ( 10) tTC ( 10) 

INTEGER  AA.BB 
REAL*8  F  I L X 
REAL*8  FILY 
REAL*8  FILR 
REAL*8  FILT 

THE  ARRAYS  MAY  BE  DIMENSIONED  TO  HANDLE  ANY  REASONABLE 
NUMBER  OF  ANGLES.  TB  AND  TC  ARE  ARRAYS  USED  FOR  PRINTING 
OUT  DATA. 

FOR  INPUT, USUALLY  LET  LA1-LB 1* I NCA* INCB=1 .  LA2*NUMBER 
OF  PITCH  ANGLES  AND  LB2-NUMBER  OF  HORIZONTAL  ANGLES  REQUESTED. 
ARESOL  AND  BRESOL  ARE  THE  RESOLUTIONS  FOR  PITCH  AND  HORIZ. 
ANGLES,  RESPECTIVELY.  ASTART  AND  BSTART  ARE  THE  STARTING  POINTS 
OF  SAME. 

NX  IS  THE  NUMBER  OF  ROWS  OF  MODULES  AND  NY  IS  THE  NUMBER 
OF  COLUMNS. 

NORML  IS  THE  FLAG  FOR  NORMALIZATION. 

ROLL  IS  THE  ROLL  ANGLE  TO  BE  PLOTTED.  VWEAP  IS  THE 
WEAPON'S  VELOCITY(IN  KNOTS),  AND  AMAINR  AND  BMAINR  DESIGNATE 
THE  MAIN  RESPONSE  AXIS(PITCH  AND  HORIZ.  ANGLES,  RESPECTIVELY). 

0PEN(UNIT  =  25,TYPE*,0LD,',NAME*'PL0TDAT') 
0PEN(UNIT=45,TYPE='0LD' , NAME- ' THR YPLT ' ) 

A0FSET=0. 

BOFSETsO. 

ANORML-O. 

DI--777. 

R0LL=0. 

P I *3. 141 59 
DTR=. 0174533 
0FFLAG=0. 

READ(45,*)PHASEPL0T 
READ (45 ,950)F ILX 

0PEN(UNIT*9,TYPEa'0LD'  ,NAME*F ILX) 

READ (45 ,950)FILY 
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0PEN(UNIT=10,TYPE='0LD',NAME=FILY) 

READ(45,950)FILR 

0PEN(UNIT  =  11,TYPE='0LDI ,NAME  =  F  ILR ) 

READ (45  , * ) TSCALC 
IF (TSCALC. EQ. 1. )  GO  TO  75 
READ (45 ,950)F I L  T 

0PEN(UNIT=12,TYPE=,0LD ' ,NAME=FILT) 

950  FORMAT (A) 

75  READ (45 ,*)LA1 ,LA2,INCA,LB1,LB2,INCB,ARES0L,BRES0L , ASTART  ,  B START 
READ(45,*)F,C,NX,NY 

FLAGS  ARE  TO  BE  ENTERED  AS  FOLLOWS:  A  #1#  MEANS  GO,  AND 
ANYTHING  ELSE  WILL  BE  TAKEN  TO  MEAN  NO  GO. 


C 

C 


C 

C 

C 

C 


READ(45,*)VWEAP,AMAINR,BMAINR,N0RML,FREQFL 

PCHEBP=AMA I  NR 

IF(N0RML.NE.2)G0  TO  70 

READ(45,*)AN0RML 

AN0RML=10**(AN0RML/20) 

70  CONTINUE 

READ (45 ,*)OFFLAG 

IF (OFFLAG.NE. 1)G0  TO  78 

AOFSET=AMAINR 


BOFSET=BMAINR 
78  CONTINUE 

READ (45 ,*)FROLL 
IF(FR0LL.NE.1)G0  TO  71 
READ (45 ,*)ROLL 

71  CONTINUE 

READ  (45  ,*)HFLAG 

READ  (9,*)((X(I,J),J=1,NX),I=1,NY) 

READ  (10,*W(Y(I,J),J  =  1,NX)  ,1  =  1, NY) 

READ  (11,*)((R(I»J).J*1.NX),I=1,NY) 

IF (TSCALC. EQ.O)GO  TO  72 

CALL  PHASTR ( X , Y ,R ,T ,NX ,N Y , AMAINR  ,BMAINR  ,P I  ,DTR  ) 

GO  TO  73 

72  CONTINUE 

READ  (12,*)((T(I,J),J=1,NX),I=1,NY) 

73  CONTINUE 

CHANGE  VWEAP  FROM  KNOTS  TO  FT/SEC 
VWEAP=VWEAP*1.689 
COMPUTE  NUMBER  OF  ENABLED  ARRAYS 
RT  =  0 

DO  5  1  =  1, NY 
DO  4  J*1  ,NX 
RT=RT+R(I,J) 

4  CONTINUE 

SUMIN1TANDESUMIN2  ARE  THE  VOLUME  AND  SURFACE (BOTTOM )  EQUIV. 
BEAM  PATTERNS.  INTEGRATIONS  ARE  DONE  BY  SIMPLY  SUMMING  UP 
OVER  THE  VOLUME  OR  SURFACE  IN  INCREMENTS  OF  THE  APPROPRIATE 
RESOLUTION(S) .  THE  VOLUME  INTEGRAL  MUST  BE  DONE  OVER  THE  E 


o  o  o  o  o 
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C  TIRE  VOLUME  TO  HAVE  ANY  MEANING,  AND  THE  SURFACE ( BOTTOM )  IN- 
C  TEGRAL  IS  TAKEN  AT  THE  LARGEST  POINT,  USUALLY  0  DEGREE  PITCH. 
SUMIN1=0 
SUM  I N2  =  0 
D  1 1  NT  =  0 

W=(2*PI*F*(C+VWEAP)  *1000)/ (CM2*  (C-VWEAP) ) 

IF  (NORML . NE . 1 )GO  TO  1 

C  CALCULATE  THE  NORMALIZATION  VALUE  ( ANORML ) 

IF (FROLL. EQ. 0)G0  TO  74 

IF (FROLL.  EQ.  1 )  AMAINR-( AS IN(S IN (ROLL*DTR )*S IN (BMAINR*DTR ) )  )/DTR 
IF (FROLL.  EQ.  1 )BMA INR  =  ( ATAN ( TAN (BMAI NR*DTR ) *COS ( ROLL*DTR  ) ))/DTR 
IF(FROLL.EQ.l.AND.OFFLAG.NE.l)GO  TO  74 
AMAINR=AOFSET 
BMAINR=BOFSET 
74  CAMR=COS ( AMAI NR*DTR ) 

SAMR=SIN (AMAINR*DTR ) 

CBMR=COS(BMAINR*DTR) 

SBMR=SIN(BMAINR*DTR) 

E-l 

H=((.6+CAMR*CBMR)**2)/2.56 
IF  (HFLAG.EQ. 1 ) H  =  1 
IF (H. LT. 0. )  H  =  -H 
RRT  =  0 
A 1 1 T  =  0 

DO  611  1  =  1, NY 
DO  612  J  =  1 ,  N X 

DDD=(W*(X(I  ,J)*SBMR*CAMR  +  Y(I  ,J)*SAMR)+T(I  ,J)*DTR) 
RR=R(I,J)*COS(DDD) 

AII=R(I,J)*SIN(DDD) 

RRT=RRT+RR 
AI IT  =  AI IT+AI I 
612  CONTINUE 
611  CONTINUE 

AN0RML=H*SQRT(RRT**2  +  AIIT**2) 

IF (ANORML. LT.O. 00001 )AN0RML=0. 00001 
GO  TO  2 

1  IF (NORML. EQ.O)ANORML=RT 
C 

2  CONTINUE 

START  THE  MAIN  CALCULATIONS.  THE  AA  LOOP  RUNS  THROUGH  THE 
PITCH  ANGLES;  THE  BB  LOOP  THROUGH  THE  HORIZ.  ANGLES.  DB  VALUES 
LOWER  THAN  -50  ARE  SET  TO  -50. 

DO  50  AA  =  LA1  ,LA2  ,  INCA 
A(AA)=ASTART+(AA-1)*ARES0L 
PH Y= ( A ( A A )  +  AOFSET)*DTR 
CA=COS (PH Y ) 

SA=S I N ( PH Y ) 

DO  60  BB  =  LB1  ,LB2 , INCB 
B(BB)=BSTART+(BB-1)*BRES0L 
TS  I ■ (B (BB )  +  BOFSET)*DTR 
SB-SIN(TSI) 


CB=COS (TS  I ) 

C  MODULE  EFFICIENCY 
C  CHECK  FOR  ROLL  ANGLE  DATA 
IF (FROLL. EQ. 0 ) G 0  TO  10 
PHY=ASIN(SIN(ROLL*DTR)  *  S IN ( B (BB ) *DTR ) ) 
TSI=ATAN(TAN(B(BB)*DTR)  *  COS (ROLL*DTR  ) ) 
IF(B(BB).LE.-90.0R.B(BB).GE.90)TSI=TSI+180*DTR 
IF(R0LL.EQ.90.0R.R0LL.EQ.-90)TSI=0 
PH Y=PHY  +  AOFSET*DTR 
TS  I  =  TS I  +  BOFSET*DTR 
SB  =  S  IN  (TS  I ) 

C8=C0S (TS I ) 

SA  =  S I N ( PH  Y ) 

CA=COS(PHY) 

10  CONTINUE 
E  =  1 

H=((.6+CA*CB)**2)/2.56 
IF (HFLAG. EQ. 1 )  H  =  1 
IF  (H.LT.O.  )H  =  -H 
RRT  =  0 
A 1 1 T=0 

C  ADJUST  FOR  FREQUENCY  ANGLE  DEPENDENCE  IF  REQUIRED 
IF (FREQFL. N E • 1 ) GO  TO  76 
VWTEMP=VWEAP*CA*CB 

W=(2*PI*F*(C+VWTEMP)*1000)/(C*12*(C-VWTEMP)) 

C  CALCULATE  AND  SUM  THE  CONTRIBUTIONS  OF  THE  INDIVIDUAL  ELEMENTS. 
C  THIS  SECTION  USES  EULER "S  THEOREM  TO  ADD  IN  THE  PHASE  SHADING. 
76  DO  61  1*1  ,NY 
DO  62  J  =  1  ,NX 

DD=(W*(X(I  ,J)*SB*CA+Y(I ,J }*SA  )+T  ( I ,  J )*DTR  ) 

RR=R ( I ,J )*COS(DD) 

A 1 1 *R ( I  , J ) *S I N (DD  ) 

RRT =RRT+RR 
AIIT*AIIT+AII 
62  CONTINUE 
61  CONTINUE 

C  TAKE  THE  MODULUS  OF  THE  COMPLEX  VALUE  (COS(X)  +  IS  IN(X ) ) 

AMOD  =  H*SQRT (RRT*RRT  +  A I IT*A I  IT) 

IF (AMOD.LT.O. 00001 )AM0D=0. 00001 
VV-20*AL0G10(AM0D/AN0RML) 

IF(VV.LT.-50)VV*-50 
IF(AI IT. EQ. 0)G0  TO  761 
ROVERI*ABS(RRT/AIIT) 

GO  TO  762 

761  ROVER  1*100. 

C  DETERMINE  THE  SIGN 

762  S*1 

IF{R0VERI.GT.0.1.AND.R0VERI.LT.10)S*2 

IF(R0VERI.LE.0.1)S-3 

IF(RRT.LT.O)S*-S 

IF (VV. LE. -50)S*0 

PHANG-ATAN2(AIIT,RRT)/DTR 

C  THIS  GIVES  THE  INTENSITY  RELATIVE  TO  THE  NORMALIZATION  FACTOR 
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A INTEN* (AMOD/ANORML )  **2 
DIINT*DIINT+CA*A  INTEN 
SUMIN1=SUMIN1+CA*AINTEN**2 

IF((PHY/DTR).£Q.PCHEBP)SUMIN2=SUMIN2+AINTEN**2 
C  WRITE  DBS  AND  SIGNS  TO  A  STORAGE  FILE  FOR  USE  IN  PRINTING 
C  AND  PLOTTING  THE  DATA  AT  A  LATER  TIME  IN  THE  PROGRAM. 

600  WR ITE ( 1 3  ,* )  VV , S  ,PHANG 
60  CONTINUE 
50  CONTINUE 

SUMIN1=SUMIN1*ARES0L*BRES0L*DTR*DTR 
SUMIN2=SUMIN2*BRES0L*DTR 
DI INT=DI INT*ARESOL*BRESOL*DTR*DTR 
IF(SUMIN2.LT..00001)SUMIN2=. 00001 
IF (SUM  INI. LT. .00001 )SUM INI*. 00001 
IF (Oil  NT. LT. .00001)0 1 1  NT*. 00001 
01  =  10.9921  -  10*AL0G 10 ( D 1 1  NT ) 

SUMIN1=10*AL0G10(SUMIN1) 

SUMIN2=10*AL0G10(SUMIN2) 

25  REWIND  13 
K 1=LB 1 
K2-LB2 
Q  =  LA2 
P  =  LA1 

IF (PHASEPL0T.NE.3)G0  TO  40 

C  READ  DATA  FROM  TAPE13  AND  WRITE  IT  OUT  TO  TAPE2.  THIS  TAPE2 
C  WILL  CONTAIN  THE  DATA  IN  AN  EASY-TO-READ  FORM. 

WRITE (2 ,110) 

163  I V= ( LB2-LB 1 ) / I NCB 

IF  (IV.LT.O.)  GO  TO  170 
IF  ( IV.LE.45. )  GO  TO  125 
LB2=LB 1+45*1 NCB 
125  L=((LA2-LA1)/INCA) 

I  L=0 

IF  (L.LT.O.)  GO  TO  160 
IF  (L.LE.4.)  GO  TO  140 
I L  *  1 

LA2«LA1+4*INCA 

140  WR  ITE  ( 2,93 )  (A(AA)-,AA  =  LA1  ,IA2 ,  INCA) 

WR I TE ( 2 , 97 ) 

WR ITE ( 2 , 97 ) 

MINC*(LA2-LA1)/INCA+1 
DO  145  M*1  ,MINC 
TB (M)*8HDB  +/- 

145  CONTINUE 

WRITE(2,109)(TB(M),M=1,MINC) 

WRITE (2,97) 

READ (13,*,  END-1 7)  ( V 1 ( BB ) ,SS 1 ( BB )  ,  XXX  ,BB  =  K  1  ,K2 , 1 NCB  ) 

RE AD (13,*, END *512)  ( V2 (BB ) , SS2(BB )  ,XXX  ,BB*K1  ,K2 , 1 NCB  ) 

GO  TO  513 

512  DO  21  BB*K1 ,K2  ,  INCB 

WRITE (2, 95)  B(BB),V1(BB),SS1(BB) 

21  CONTINUE 
GO  TO  17 


E-5 


o  o  o 


513  READ(13,*,  END-514)  ( V3 (BB  )  ,SS3(BB ), XXX ,BB=K1 ,K2 , INCB ) 

GO  TO  515 

514  DO  22  BB*K1  ,K2 , INCB 

WRITE (2, 95)  B(BB) , VI (BB ) , SSI ( BB  ) ,V2(BB) ,SS2(BB) 

22  CONTINUE 
GO  TO  17 

515  READ(13  ,*,END«516)  ( V4(BB )  ,SS4 (BB ), XXX ,BB*K1 ,K2 , INCB ) 

GO  TO  517 

516  DO  23  BB*K1,K2,INCB 

WRITE (2, 95)  B(BB)  , Vl(BB)  ,SS1(BB)  ,V2(BB) ,SS2(BB) ,V3(BB) ,SS3(BB) 

23  CONTINUE 
GO  TO  17 

517  READ(13,*,END  =  518)  ( V5 ( BB  )  , SS5 ( BB )  , XXX ,BB *K1 ,K2 , INCB ) 

GO  TO  519 

518  DO  24  BB  =  K1  ,K2  ,  INCB 

WRITE  (2, 95)  B(BB)-tVl(BB),SSl(BB)  ,V2(BB)  ,SS2(BB)  ,V3(BB)  ,SS3(BB) 
*V4(BB) ,SS4(BB) 

24  CONTINUE 
GO  TO  17 

519  DO  16  BB*K1 ,K2  ,  INCB 

WRITE (2, 95)  B ( BB )  , V 1 ( B B )  ,SS1(BB)  ,V2(BB) ,SS2(BB) ,V3(BB) ,SS3(BB) 
* V  4 ( B  B ) ,SS4(BB),V5(BB) , SS5 ( BB ) 

16  CONTINUE  , 

17  CONTINUE 
LA1=LA2+INCA 
LA2-Q 

IF  (LA1.GT.LA2)G0  TO  170 
IF  (IL.EQ.l.)  WR ITE (2 , 1 10) 

GO  TO  125 
160  LA1=P 
165  WR ITE ( 2 , 1 10 ) 

LB1«LB2+INCB 
LB2*K2 
GO  TO  163 
26  REWIND  13 
170  LB1=K1 
LB2=K2 
LA1-P 
LA2-Q 

40  CONTINUE 
C 

C  CHANGE  VWEAP  BACK  TO  KNOTS 
VWEAP=VWEAP/1.689 

WRITE  INPUT  DATA  TO  TAPE3 

AN0RML»20*AL0G10( ANORML ) 

WRITE (3,110) 

WRITE (3, 98) (LAI, LA2, INCA, ARESOL , ASTART ) 
WRITE(3,99)(K1,K2,INCB,BRES0L,BSTART) 

WRITE (3, 91)FR0LL, ROLL 
WRITE(3,114)(F,C,NX,NY) 

WR ITE ( 3 , 104  JAM A I  NR ,BMA INR .VWEAP 
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WRITE (3,1 23 )PHASEPLOT,FREQFL 
WRITE(3,106)N0RML,AN0RML 

106  FORMAT ( 1 HO  , ' NORML-  1 1 , 5X  NORMALIZATION  FACTOR  (DB)-',F10. 

WRITE(3, 107)01 

107  F0RMAT(1H0, 'DIRECTIVITY  INDEX- '  ,F10. 4) 
WRITE(3,108)A0FSET,B0FSET,HFLAG,TSCALC 

108  FORMAT ( 1H0 , 1 A  OFFSET- '  ,F8. 2 ,5X  B  OFFSET- *  ,F8. 2 ,5X  ,' HFLAG- '  , 
*F5.0,5X,'TSCALC»' ,F5.0) 

WRITE (3, 105) SUM  INI  .SUMIN2 

104  FORMAT (1 HO ,  'AMAINR-  ' , F 7 . 2 , 3X  , '  BMA I  NR-  ',F7.2,3X, 

*  1  VWEAP  (KNOTS)-  ',F7.2) 

105  FORMAT ( 1H0 , 1 EQ  BEAM  PATTERN-VOL- ' , F10. 4 , 5X  , 

* 1 EQ  BEAM  PATTERN-S.B-' ,F10.4) 

622  FORMAT ( 1H  ) 

623  FORMAT ( 1H  , 9( F 10. 3 , 3X ) ) 

45  F0RMAT(1X,6I6,4(2X,F7.2)) 

91  F0RMAT(1H0, 'ROLL  FLAG-  ' ,F3. 0 , 5X ,' ROLL  ANGLE- '  ,F8. 2  , 

*2X,  'DEGREES' ) 

93  FORMAT ( IX  , '  P ITCH  ' , 4X ,8 ( F6. 2 , 9X ) ) 

95  FORMAT ( 1H  ,5X ,F6. 2 , 3X ,8 (F5. 1 ,3X ,F3. 0 ,4X ) ) 

97  FORMAT ( 1H  ) 

123  F0RMAT(1H0, 'PLOTTING  FLAG-  '  ,F3. 0 , 5X  , ' FREQFL-  \F3.0) 

98  F0RMAT(1X,'LA1»' ,I4,4X,,LA2«',I4,4X,'INCA»' ,I4,4X,'ARES0L«' , 
*F5. 2 , 3X , ' ASTART- ‘  ,F7. 2 ) 

99  F0RMAT(1X,'LB1-' ,I4,4X,'LB2«' ,I4,4X,' INCB-' ,I4,4X,'BRES0L-' , 
*F5.2,3X,'BSTART«' ,F7.2) 

109  FORMAT ( IX ,2X , '  AZ  \4X,8(2X,A8,5X)) 

110  FORMAT(lHl) 

114  FORMAT (1H0,1 FREQ. ( KHZ ) - ' ,F8. 3 , 5X  , ' SND.  SPD. (F/S ) - ' , F5. 0 , 5X , 
l'NX-' ,I4,5X, 'NY-'  ,  14 ,5X) 
IF(PHASEPL0T.EQ.0.0R.PHASEPL0T.E0.1)CALL  PLOTDB 
IF(PHASEPL0T.EQ.1.0R.PHASEPL0T.EQ.2)CALL  PLOTPHANG 
CLOSE (UN IT-9, DISPOSE -'SAVE' ) 

CLOSE (UN IT-10, DISPOSE-' SAVE' ) 

CLOSE (UN IT-11, DISPOSE -'SAVE  ' ) 

CLOSE (UN IT-12, DISPOSE-' SAVE'  ) 

CLOSE (UN IT-25, DISPOSE -'SAVE' ) 

STOP 

END 


i 


c 

c 


SUBROUTINE  PHASTR ( X , Y , R ,T  . NX  , N Y  , AMA INR  ,BMA INR  ,P I ,DTR ) 
DIMENSION  X(8,8),Y(8,8),R(8,8),T(8,8) 

READ(45,*)FD,CD 
SA-S IN(AMA I NR*DTR ) 

SB«SIN(BMAINR*DTR ) 

CA»COS(AMAINR*DTR ) 

CB»COS(BMAINR*DTR ) 

W-(2*PI*FD*1000)/(CD*12) 

DO  10  I  - 1 ,  N  Y 
DO  20  J-l  ,NX 

T(I,J)--(W*(X(I,J) *CA*SB 
IF(R(I,0).EQ.0)T{I,J)-0 


Y ( I , J ) *SA ) ) /DTR 
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20  CONTINUE 
10  CONTINUE 
DO  30  1*1 , NX 

WRITE(20,50)(T(I,J),J«1,NY) 
50  F0RMAT(1X,9(F8.3,1X)) 

30  CONTINUE 
RETURN 
END 


SUBROUTINE  PLOTS (VI  ,SS1 ) 

DIMENSION  VI (365) ,SS1 (365) 
C0MM0N/DBSIGN/DBDAT(51,181)  ,SDAT(3,181) 
DO  10  N  =  1 ,181 
I«IABS(INT(V1(N)-.5J)  +  1 
S  =  0 

IF(SS1(N).GT.0)S=1 
IF(SSl(N).LT.O)S=-l 
K  =  S  +  2 

DBDAT ( I ,N )=1H* 

SDAT (K,N)=1H* 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  PLOTDB 

COMMON/DBS IGN/DBDAT (51,181),SDAT(3,181),  Vl(365)  ,SS1 ( 365 )  ,B ( 365 ) 
DO  680  J * 1, 1 8 1 
DO  681  1=1,51 
DBDAT ( I ,  J  )  =  1H 

681  CONTINUE 

DO  682  K-1,3 
SDAT (K  ,  J  )  =  1H 

682  CONTINUE 
680  CONTINUE 

REWIND  13 
DO  1  1=1,181 
B ( I ) -1-91 
1  CONTINUE 
WR ITE (2 ,377  ) 

WRITE (2 ,622) 

377  F0RMAT(1H1 ,8X,5( '  ANGLE  DB  SIGN* , 5X ) ) 

622  FORMAT ( 1H  ) 

READ(13,*-,END*515)(V1(BB),SS1(BB)-,XX,BB  =  1,181) 

515  CONTINUE 

CALL  PLOTS (VI  ,SS1 ) 

I MM* 181/5 
MM-IMM+1 
DO  517  1  =  1, MM 

WRITE (2,516)  (B(5*(I-1)+J)-,V1(5*(I-1)+J)-, SSI (5*(I-1)+J)  ,J  =  1,5) 
517  CONTINUE 

516  F0RMAT(1H  ,8X,5(F6.1,3X,F5.1,2X,F3.0,5X)) 
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26  REWIND  13 
WRITE (25, 650) 

650  F0RMAT(1H1,87X, 'DECIBELS  VS  ANGLE') 
WR ITE ( 25 ,622) 

WRITE (25, 651)(DBDAT(1,J),J  =  1, 181) 

DO  671  1*2,10 

WRITE (25, 652) (DBDAT ( I ,0)  ,0  =  1,181) 

671  CONTINUE 

651  FORMAT  ( 1H  ,2X,'  0  '  ,2X  181  (A1 ) ) 

652  FORMAT  ( 1H  , 7X , '  +  '  ,  18 1 ( Al ) ) 

WRITE (25, 653) (DBDAT (11,0) ,0=1,181) 

653  FORMAT ( 1H  ,2X , ' -10 ' , 2X , '  +  '  ,181 (A1 ) ) 

DO  672  1*12,20 

WR ITE (25,652) (DBDAT ( I  ,0)  ,0  =  1, 181) 

672  CONTINUE 

WR ITE ( 25, 654 ) (DBDAT ( 2 1,0),0* 1,181) 

654  FORMAT ( 1H  , 2X , ' -20 ' , 2X  , '  +  ' , 18 1 ( A1 ) ) 

DO  673  1*22,30 

WRITE(25,652)(DBDAT(I ,0 ),0*1 , 181 ) 

673  CONTINUE 

WRITE(25,655)(DBDAT(31,0),0=1,181) 

655  FORMAT ( 1H  ,2X, '-30' ,2X, '+' ,181 (Al) ) 

DO  674  1*32,40 

WRITE  (25, 652)  (DBDAT  (I ,  J  ) -,  J*  1 , 181 ) 

674  CONTINUE 

WRITE  (25, 656)  (DBDAT (41, 0),0=1, 181) 

656  FORMAT ( 1H  , 2X , ' -40 ' ,2X , ' + ' , 181 ( Al ) ) 

DO  675  1*42,50 

WRITE (25, 652) (DBDAT (I  ,0)  ,0=1,181) 

675  CONTINUE 


WR ITE ( 25, 657 ) (DBDAT ( 51, J),  0  =  1, 181) 

657  FORMAT ( 1H  ,2X , ' -50 '  ,2X  , '  +  '  ,181 ( Al ) ) 

WR ITE (25 ,658) 

658  FORMAT ( 1H  ,8X,18('+ . ')>,'  +  ') 

WR ITE (25 ,659) 

659  FORMAT (1H  ,6X , ' -90 ' , 7X , ' -80 ' , 7X , ' - 70 '  , 7X , ' -60 '  , 7X , ' -50 '  , 
*'-40'  ,7X,'-30‘  ,7X  , ' -20 '  , 7 X  ,  *  - 1 0 '  ,9X  ,'0' ,8X,'10' ,8X,'20' , 
*8X , '  30 '  ‘,8X't '  40 '  -,8X', '  50*  >,8X  , ' 60'  <,8X , '  70*  ',8X , ' 80 '  -,8X , '  90 '  ) 

WR ITE (25 ,679) 

679  FORMAT ( 1H  ,/////) 

WR ITE ( 25 ,660) 

660  FORMAT ( 1H  ,89X,'SIGN  VS  ANGLE') 

WRITE (25,622) 

WRITE (25, 661) (SDAT (3,J),J=1,181) 

661  FORMAT ( 1H  ,3X,‘  1 ' , 2X , ' + ’ ,181 ( Al ) ) 

DO  676  1-1,7 

WRITE(25,662) 

676  CONTINUE 


662  FORMAT ( 1H  , 7 X , '  +  ' ) 

WRITE (25, 663) (SDAT (2,0) ,0-1 ,181) 
DO  677  1-1,7 
WR ITE (25,662) 

677  CONTINUE 


7X , 


E-9 


NSWC  TR  82-376 


663  FORMAT ( 1H  ,3X,‘  0 1 ,2X , *  +  * ,181 ( A1 ) ) 

WRITE (25, 664) (SDAT (1,J),J®1,181) 

664  FORMAT ( 1H  ,  3X  , 1 -  1 ' , 2X  ,  •  +  1  , 1 81 ( A1 ) ) 
WRITE(25,658) 

WR ITE (25 ,659) 

END 

C 

SUBROUTINE  PLOTPHANG 

DIMENSION  SSI (365)  .DBPHASE ( 73 , 181 ) , B ( 3 6 5 ) 
REWIND  13  ' 

RE  AD  (13  #*,END*515)(X,Y,SS1(I)-,  1  =  1,181) 

515  CONTINUE 
DO  1  1=1,181 
B(  I  )  =  I-91 

1  CONTINUE 

DO  2  J  =  1 ,181 
DO  2  1=1,73 
DBPHASE ( 1 ,0 )  =  1H 
2  CONTINUE 

781  DO  10  N  =  1 ,181 

I=INT( (-SS1 (N)/5. )+37. 5) 

IF  ( I . GT. 73)1  =  73 
IF(I.LT. 1)1=1 

DBPHASE ( I ,N ) = 1H*  ' 

10  CONTINUE 

C  WRITE  OUT  NICE  EASILY  READ  DATA  SET 
WRITE (2, 377) 

377  FORMAT  (1H1 ,8X  ,5( '  ANGLE  PHASE'. ,9X)) 

IMM=181/5  7 

MM-IMM+1 

DO  517  1=1, MM 

WRITE(2,516)(B(5*(I-1)+0),SS1(5*(I-1)+0),J=1,5) 

516  FORMAT ( 1H  ,8X, 5(F6. 1 ,3X ,F6. 1 ,9X) ) 

517  CONTINUE 

WRITE (25,650) 

650  F0RMAT(1H1 ,87X, ‘PHASE  ANGLE  VS  THETA') 

WR ITE (25 ,622) 

622  FORMA! (1H  ) 

WR ITE (25, 651 ) (DBPHASE (1,J),J= 1,181) 

DO  671  1-2,9 

WR  ITE  (25, 652)  (DBPHASE  ( I,  J).,J.  1,181) 

671  CONTINUE 

651  FORMAT ( 1H  ,2X , ' 1 80 ' , 2X , ' + » , 181 ( A1 ) ) 

652  FORMAT ( 1H  ,7X  , ’  + '  ,181  (A1 ) ) 

WR ITE (25, 653) (DBPHASE ( 10,  J),J  =  1, 181) 

653  FORMAT ( 1H  , 2X , ' 1 35 ' , 2X , ' + • , 181 ( A1 ) ) 

DO  672  1=11,18 

WR ITE (25, 652) (DBPHASE ( I, J),J. 1,181) 

672  CONTINUE 

WR  ITE  (25, 654)  (DBPHASE  ( 19,  J).,J  =  1, 181) 

654  FORMAT (1H  ,2X,‘  90' ,2X, '  +  '  ,181(A1)) 

DO  673  1=20,27 

WR ITE (25, 652) (DBPHASE ( I, J),J  =  1, 181) 
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673  CONTINUE 

WR ITE ( 25,655) (DBPHASE ( 28 , J )  ,J  =  1,181) 

655  FORMAT  ( 1H  ,2X,'  45 ' ,2X , '  +  '  ,181 (A1 ) ) 

DO  674  1=29,36 

WRITE (25, 652) (DBPHASE (I, J)  ,J  =  1,181) 

674  CONTINUE 

WR ITE (25, 656) (DBPHASE ( 37, J)  ,J*1,181) 

656  F0RMAT(1H  ,2X,'  0 ' ,2X , * + ' , 181 ( A1 ) ) 

DO  675  1=38,45 

WRITE (25,652) (DBPHASE ( I ,J),J=1,181) 

675  CONTINUE 

WRITE (25, 657) (DBPHASE (46,J),J  =  1, 181) 

657  FORMAT ( 1H  , 2X , ' -45  * ,2X , '  +  ‘ , 181 ( A1 ) ) 

DO  845  1=47,54 

WR ITE ( 25, 652 ) (DBPHASE ( I, J),J  = 1,181) 

845  CONTINUE 

WRITE (25, 847) (DBPHASE (55, J)  ,J  =  1, 181) 

847  FORMAT ( 1H  ,2X , ' - 90 ' , 2X , *  +  ' ,181 ( A1 ) ) 

DO  848  1=56,63 

WRITE (25, 652) (DBPHASE (I ,J),J= 1,181) 

848  CONTINUE 

WR ITE ( 25, 849) (DBPHASE (64,J),J= 1,181) 

849  FORMAT  ( 1H  ,  1X-, '  - 1 35 '  ,2X  ,  *  +  '  ,  181  (  A  1 ) ) 

DO  850  1=65,72 

WRITE (25, 652) (DBPHASE (I , J )  ,J  =  1,181) 

850  CONTINUE 

WR ITE (25, 851 ) (DBPHASE (7?,  J),J  =  1, 181) 

851  FORMAT ( 1H  , 1 X , ' -1 80 ' ,2X , '  +  ' , 181 ( A1 ) ) 

WRITE (25 ,658) 

658  FORMAT ( 1H  ,8X,18('+ . '),'+') 

WR ITE (25 ,659) 

659  FORMAT ( 1H  ,6X  , ' -90 1  ,7X , ' -80  ' -,7X , ' -70  '  ,7X , ' -60  ' 
*'-40' ,7X,'-30' , 7X , ' - 20 ' ,7X, '-10' ,9X, '0' ,8X,  '10 
*8X  , 1  30  1  .eX.MO'/.SX.'SO-  ,8X,  '60*  .SX.^O'  ,8X,'80 

1000  FORMAT ( 1H  ,13) 

RETURN 
END 
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